Plots

#install.packages("BiocManager")
#BiocManager::install("Seurat")
#BiocManager::install("ggplot2")
#BiocManager::install("sctransform")

library(Seurat)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ape)
library(cowplot)
library(Matrix)
library(EnhancedVolcano)

## Loading required package: ggrepel

#find out where you are
getwd()

## [1] "/Users/whippoorwill/Desktop/Code/scrnaseq/Rmd"

#Specify where your matrix files are
dir= "/Users/whippoorwill/Desktop/Sequencing/LD_AVM02/"
datafolder = "Data/Seurat"
filename = "Microglia_BC_Macrophages_subset.RData"
organism = "Mouse"
defile = "Macrophage_only_all_markers.csv"

m = c("nCount_RNA","nFeature_RNA","percent.mito","percent.ribo")
Plotfolder = "Plots"

if(organism == "Mouse"){library(org.Mm.eg.db)}

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following object is masked from 'package:Matrix':
## 
##     which

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:Matrix':
## 
##     expand

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:base':
## 
##     expand.grid

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:dplyr':
## 
##     select

## 

if(organism == "Human"){library(org.Hs.eg.db)}
if(organism == "Zebrafish"){library(org.Dr.eg.db)}

Load in your filtered dataset

load(file.path(dir,datafolder,filename))

Other differential expression methods:

1. Compare each cluster to one “baseline” cluster

base = '0'                 #default cluster
column = "seurat_clusters" #can pick any metadata column

Idents(sobject) = column
clusters = levels(sobject@meta.data[,column]) 
clusters = clusters[clusters != base]

for (cluster in clusters){
  markers_all <- FindMarkers(
    object = sobject,
    ident.1 = cluster,
    ident.2 = base,
    only.pos = FALSE, 
    min.pct = 0.10, #gene must be present in 10% of the cells in the cluster
    logfc.threshold = 0,
    test.use = "MAST")
  dim(markers_all)
  head(markers_all)
  write.csv(markers_all,file = file.path(dir,"Spreadsheets",paste0("markers",cluster,"_vs_",base,".csv")))
}

Or you can run “bulk-seq” like analyses based on your original sample IDs:

#can pick any metadata column
column = "sample_description"
#default cluster
cluster1 = 'Deprived-P5'    
#cluster of interest
cluster2 = "Control-P5"

Idents(sobject) = column
markers = FindMarkers(sobject,
                      ident.1=cluster1,
                      ident.2 = cluster2,
                      only.pos=F,
                      logfc.threshold = 0.0,
                      min.pct = 0.1,
                      test.use = "MAST")
write.csv(markers,file.path(dir,"Spreadsheets",paste0("markers_",cluster1,"_vs_",cluster2,".csv")))

You can also make plots with only subsets of cells: Best for making umap plots where both conditions/individuals have equal representation.

#first check how many cells are in each group so you don't pick a number more than the min
column = "sample_description"
ncells = 2000
genes = c("Cx3cr1","P2ry12","Spp1","Nrxn2","Ifitm3","Mki67","Pf4")
features = m
categories = c("sample_description","age","condition","celltypecluster")
name = "macrophage_equalcells"

Function to print multiple graphs:

PrintSeuratGraph = function(namecard = "a",seurat_object = sobject,graphtype = "feature",feature = NULL,group = NULL,split=NULL,cellnames=NULL){
  if (!is.null(cellnames)){
    Idents(seurat_object) = cellnames[1]
    cells = colnames(seurat_object)[Idents(seurat_object) %in% cellnames[2:length(cellnames)]]} 
  else {cells = cellnames}
  if (graphtype == "feature"){
    graph = FeaturePlot(seurat_object,features = feature,split.by = split, cells = cells,cols = c("lightyellow","darkred"))
  }
  if (graphtype == "violin"){
    graph = VlnPlot(seurat_object,features = feature, pt.size = 0.1, idents = cellnames[2:length(cellnames)],group.by = group, split.by = split)
  }
  if (graphtype == "dim"){
    graph = DimPlot(seurat_object,cells = cells, group.by = group, split.by = split)
    
  }
  name = paste0(feature,"_",graphtype,namecard,".eps")
  graph
  setEPS()
  postscript(file.path(dir,Plotfolder,name))
  print(graph)
  dev.off()
}

table(sobject$sample_description)

## 
##  Control-P5  Control-P7 Deprived-P5 Deprived-P7 
##        2975        4125        2566        2743

cellnames = sobject@meta.data[,column]
names(cellnames) = colnames(sobject)
groups = levels(cellnames)

newcellnames = NULL
for (group in groups){
  cells = sample((cellnames)[cellnames == group],ncells)
  newcellnames = c(newcellnames,cells)
}

#feature plots
for(feature in c(genes,features)){
  PrintSeuratGraph(namecard = name,graphtype = "feature",feature = feature,cellnames = newcellnames)
}

#split feature plots by individual
for(feature in c(features)){
  PrintSeuratGraph(namecard = paste0(name,"_split"),graphtype = "feature",feature = feature,split = "sample_description",cellnames = newcellnames)
}

#dim plots for clustering
for(group in categories){
  PrintSeuratGraph(namecard = name,graphtype = "dim",group = group, feature = group,cellnames = newcellnames)
}

#violin plots
for(feature in c(genes,features)){
  PrintSeuratGraph(namecard = name,graphtype = "violin",feature = feature,group = "seurat_clusters",cellnames = newcellnames)
}

You can save umap and pc embeddings to load into other programs like Loupe and anything in Python like scvelo (velocity)

#Edit name only if one project has multiple embeddings
name = paste0(Project(sobject),"_",name)

#Save umap embeddings
umap =as.data.frame(sobject@reductions$umap@cell.embeddings)
umap$Cellname = rownames(umap)
umap = umap[,order(colnames(umap))]
write.csv(umap,file = file.path(dir,"Annotation",paste0(name,"_umap_embed.csv")),row.names = F)

#save pcs
pc =as.data.frame(sobject@reductions$pc@cell.embeddings)
pc$Cellname = rownames(pc)
pc = pc[,order(colnames(pc))]
write.csv(pc,file = file.path(dir,"Annotation",paste0(name,"_pc_embed.csv")),row.names = F)

#save metadata
meta = sobject@meta.data
meta$Cellname = rownames(meta)
meta = meta[,c(ncol(meta),1:(ncol(meta)-1))]
write.csv(pc,file = file.path(dir,"Annotation",paste0(name,"_metadata.csv")),row.names = F)

#For PanoView (in Python, a way to mathematically validate # of clusters found)
annotation = sobject$sample_description #or any column
x = as.matrix(GetAssayData(sobject,slot = "counts"))
x = x[rownames(x) %in% VariableFeatures(sobject),]

write.csv(x,file = file.path(dir,"Annotation",paste0(name,"_counts.csv")))
write.csv(annotation,file = file.path(dir,"Annotation",paste0(name,"_annotation.csv")))

You can add annotations from the annotationDB database:

de = read.csv(file.path(dir,"Spreadsheets",defile),stringsAsFactors = F) #any spreadsheet with gene symbols or other identifiers

if (organism == "Mouse"){db = org.Mm.eg.db}
if (organism == "Human"){db = org.Hs.eg.db}
if (organism == "Zebrafish"){db = org.Dr.eg.db}

ids=de$gene 
fromKey="SYMBOL" #must match the ids - could also be ensembl ID
toKey=c("GENENAME","ENSEMBL","UNIPROT") #whatever annotation you want to add - find with keytypes(db)
selRes<-AnnotationDbi::select(db,keys=ids,keytype=fromKey,columns=c(fromKey,toKey))

## 'select()' returned many:many mapping between keys and columns

x=selRes[match(ids,selRes[,1]),1:(length(toKey)+1)]
identical(x$SYMBOL,de$gene)

## [1] TRUE

de$GeneName = x$GENENAME
de$Ensembl = x$ENSEMBL
de$Uniprot = x$UNIPROT

Volcano Plot

Set your parameters

#Minimum fold change (i.e. 1.15 = 15% increase)
minfc = 1.15
#Max adj. p value
alpha = 1e-25
#Clusters selected
categories = levels(as.factor(de$cluster))
#Genes to highlight
ngenes = 20

Set up the spreadsheet correctly

colnames(de)[8] = "Gene"
newlist = list()
clusters = levels(as.factor(de$cluster))

#Split by cluster
i = 1
for (cluster in clusters){
  newlist[[cluster]] = de[de$cluster == cluster,]
  i = i+1
}

#select a single cluster

for (category in categories){
  fc = newlist[[category]]
  fc = fc[!is.na(fc$avg_logFC),]
  colorkeysdown = fc$Gene[fc$avg_logFC < -log2(minfc) & fc$p_val_adj < alpha]
  colorkeysup = fc$Gene[fc$avg_logFC > log2(minfc) & fc$p_val_adj < alpha]

#Either highlight specific genes or pick the top genes in colorkeysup/down
  top = fc[fc$p_val_adj<alpha,]
  top = top[order(top$avg_logFC),"Gene"]
  highlight = c(head(top,ngenes),tail(top,ngenes))

  allcolors = rep("darkgrey",length(fc$Gene))
  names(allcolors) = fc$Gene

  allcolors[names(allcolors) %in% colorkeysdown] = "blue"
  allcolors[names(allcolors) %in% colorkeysup]= "red"
  allcolors[names(allcolors) %in% highlight]= "yellow"

  names(allcolors)[allcolors == "yellow"] = "labelled"
  names(allcolors)[allcolors == "red"] = "u"
  names(allcolors)[allcolors == "darkgrey"] = "-"
  names(allcolors)[allcolors == "blue"] = "d"
  
  setEPS()
  postscript(file.path(dir,"Plots",paste0("Volcano_",category,".eps")))
  print(EnhancedVolcano(fc,
                lab = fc$Gene,
                x = 'avg_logFC',
                y = 'p_val_adj',
                xlim = c(-3, 3),
                title = category,
                subtitle = "",
                drawConnectors = F,
                legendPosition = 'right',
                legendVisible = F,
                pCutoff = alpha,
                FCcutoff = log2(minfc),
                selectLab = highlight,
                transcriptPointSize = 1.5,
                transcriptLabSize = 2.0,
                col=c('black', 'black', 'black', 'red3'),
                colCustom = allcolors,
                gridlines.major = F,
                gridlines.minor = F,
                colAlpha = 1))
  dev.off()

  print(EnhancedVolcano(fc,
                lab = fc$Gene,
                x = 'avg_logFC',
                y = 'p_val_adj',
                xlim = c(-3, 3),
                title = category,
                subtitle = "",
                drawConnectors = F,
                legendPosition = 'right',
                legendVisible = F,
                pCutoff = alpha,
                FCcutoff = log2(minfc),
                selectLab = highlight,
                transcriptPointSize = 1.5,
                transcriptLabSize = 2.0,
                col=c('black', 'black', 'black', 'red3'),
                colCustom = allcolors,
                gridlines.major = F,
                gridlines.minor = F,
                colAlpha = 1))

}

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

## Warning: One or more P values is 0. Converting to minimum possible value...

You can find what phase of the cell cycle each cell is in:

Cell cycle regression specifically for mouse genes right now - can adjust slightly for human/zebrafish data by changing the capitalization of the spreadsheet

cc.genes <- readLines(con = "../Spreadsheets/regev_lab_cell_cycle_genes.txt")
genesofinterest = c("Pcna","Top2a","Mcm6","Mki67")

#for zebrafish and mice, change gene names to lowercase
if (organism %in% c("Zebrafish","Mouse")){cc.genes<-tolower(cc.genes)}

#for mice, capitalize the first letter of each gene
if (organism == "Mouse"){
  cc.genes<-unname(sapply(cc.genes,function(x){
    x<-paste0(toupper(substr(x,start = 1,stop = 1)),substring(x,first=2))
    }))
}

#assign each gene to "S" phase or "G2M" phase
s.genes <- cc.genes[1:43]
g2m.genes <- cc.genes[44:97]

#score each cell by gene expression on this subset of genes
sobject <- CellCycleScoring(object = sobject, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

## Warning: The following features are not present in the object: Mlf1ip, not
## searching for symbol synonyms

## Warning: The following features are not present in the object: Fam64a, Hn1, not
## searching for symbol synonyms

# view cell cycle scores and phase assignments
head(x = sobject@meta.data)

##                        orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGTCTACCA-LD5LD      LD5LD      22762         4959     2.455847
## AAACGAAAGCGTGTCC-LD5LD      LD5LD      10098         3218     3.069915
## AAACGAAAGGTCTTTG-LD5LD      LD5LD      13770         4173     2.222222
## AAACGAATCGGTAGAG-LD5LD      LD5LD      19687         4451     5.338548
## AAACGCTAGGTACATA-LD5LD      LD5LD      13015         3875     3.234729
## AAACGCTAGTAGCATA-LD5LD      LD5LD      12825         3816     4.421053
##                        percent.ribo sample_description age condition    sex
## AAACCCAAGTCTACCA-LD5LD    19.299710        Deprived-P5  P5  Deprived Female
## AAACGAAAGCGTGTCC-LD5LD    22.291543        Deprived-P5  P5  Deprived   Male
## AAACGAAAGGTCTTTG-LD5LD    12.541757        Deprived-P5  P5  Deprived Female
## AAACGAATCGGTAGAG-LD5LD    20.150353        Deprived-P5  P5  Deprived Female
## AAACGCTAGGTACATA-LD5LD     9.880907        Deprived-P5  P5  Deprived Female
## AAACGCTAGTAGCATA-LD5LD    19.555556        Deprived-P5  P5  Deprived Female
##                        nCount_SCT nFeature_SCT SCT_snn_res.0.5 SCT_snn_res.1.5
## AAACCCAAGTCTACCA-LD5LD      15123         4774               4              10
## AAACGAAAGCGTGTCC-LD5LD      13442         3223               0               0
## AAACGAAAGGTCTTTG-LD5LD      13974         4172               3               2
## AAACGAATCGGTAGAG-LD5LD      15110         4409               4              10
## AAACGCTAGGTACATA-LD5LD      13714         3874               2               4
## AAACGCTAGTAGCATA-LD5LD      13672         3816               0               0
##                        SCT_snn_res.1 seurat_clusters  celltype celltypecluster
## AAACCCAAGTCTACCA-LD5LD            10               0 Microglia     Microglia-0
## AAACGAAAGCGTGTCC-LD5LD             0               0 Microglia     Microglia-0
## AAACGAAAGGTCTTTG-LD5LD             1               0 Microglia     Microglia-0
## AAACGAATCGGTAGAG-LD5LD            10               0 Microglia     Microglia-0
## AAACGCTAGGTACATA-LD5LD             8               2 Microglia     Microglia-2
## AAACGCTAGTAGCATA-LD5LD             0               0 Microglia     Microglia-0
##                        SCT_snn_res.0.1    S.Score  G2M.Score Phase old.ident
## AAACCCAAGTCTACCA-LD5LD               0 -0.3239106 -0.4315991    G1        G1
## AAACGAAAGCGTGTCC-LD5LD               0 -0.2204614 -0.3003426    G1        G1
## AAACGAAAGGTCTTTG-LD5LD               0  0.3200799 -0.3523386     S         S
## AAACGAATCGGTAGAG-LD5LD               0 -0.2900263 -0.4606724    G1        G1
## AAACGCTAGGTACATA-LD5LD               2 -0.1918543 -0.3783492    G1        G1
## AAACGCTAGTAGCATA-LD5LD               0 -0.1110625 -0.3937208    G1        G1

RidgePlot(object = sobject, features = genesofinterest)

## Picking joint bandwidth of 0.0685

## Picking joint bandwidth of 0.123

## Picking joint bandwidth of 0.0844

## Picking joint bandwidth of 0.135

FeaturePlot(sobject,"S.Score")

FeaturePlot(sobject,"G2M.Score")

#save phase plot
setEPS()
postscript(file.path(dir,"Plots",paste0(name,"_phasechart.eps")))
DimPlot(sobject,group.by = "Phase")
dev.off()

## quartz_off_screen 
##                 2

#save seurat object with cell cycle scoring
save(sobject,file = file.path(dir,datafolder,filename)))

You can regress out the phase genes

regress = c("S.Score","G2M.Score","nCount_RNA","percent.mito")
columns = c("sample_description","seurat_clusters","celltypecluster","condition","age","Phase")

#run PCA on only the cell cycle genes
sobject2 <- RunPCA(object = sobject, pc.genes = c(s.genes, g2m.genes), do.print = FALSE, maxit=10000)

#cell cycle only
PCAPlot(object = sobject2)

#original
PCAPlot(object = sobject)

#regress out cell cycle variables
sobject2 <- ScaleData(object = sobject2, features = VariableFeatures(sobject), vars.to.regress = regress)

# Now, a PCA on the variable genes no longer returns components associated
# with cell cycle
sobject2 <- RunPCA(object = sobject2, features = VariableFeatures(sobject), genes.print = 10)
PCAPlot(object = sobject2)

#remake umap with new calculations without cell cycle
sobject2@reductions$UMAP<-NULL
sobject2<-RunUMAP(sobject2,reduction = "pca",dims = 1:30, verbose = F)
sobject2<-FindNeighbors(sobject2,dims=1:30,verbose=F)
set.seed(1)
sobject2<-FindClusters(sobject2,verbose=F,resolution = 0.5)

load(file.path(dir,datafolder, paste0(name,"_cellcycleregressed.RData")))

#make plots without cell cycle
for (column in columns){
  print(DimPlot(object = sobject2, group.by=column, pt.size=0.5,label = T))
}

## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.

#save seurat object with cell cycle scoring
save(sobject2,file = file.path(dir,datafolder,paste0(name,"_cellcycleregressed.RData")))

Barplot for any two (or more) categories

#Pick metadata columns
clustercolumn = "seurat_clusters"
samplecolumn = "sample_description"
#pick a reasonable number of cells per sample to normalize by
ncells = 2000 
cols = c("blue","red")
#If you want to only compare particular samples/conditions, split further by another metadata column: 
split = "age"

splitby = levels(as.factor(sobject[[split]][,1]))

#Make a table and normalize
r = table(sobject[[clustercolumn]][,1],sobject[[samplecolumn]][,1])

#Split the table
for (i in splitby){
  x = grep(i,colnames(r))
  t = r[,x]
  
  #remove any clusters that don't have cells
  t = t[rowSums(t)>0,]
  
  #normalize by sample
  t = apply(t,MARGIN = 2,function(x)x/sum(x))
  t = round(t*ncells,0)
  
  #convert to percents for each cluster
  t = apply(t,MARGIN = 1,function(x)x/sum(x))
  t = round(t*100,2)
  
  setEPS()
  postscript(file.path(dir,"Plots",paste0(name,i,"barplot.eps")))
  barplot(t, main="Cluster composition by percent of celltype",
        xlab="Cluster", ylab = "% of cluster", ylim = c(0,100), col=cols,axisnames = T,
        width = .2,xlim = c(0,5),legend = rownames(t), space = 0.6,cex.names = 0.8,axis.lty = 1)
  dev.off()
  print(barplot(t, main="Cluster composition by percent of celltype",
        xlab="Cluster", ylab = "% of cluster", ylim = c(0,100), col=cols,axisnames = T,
        width = .2,xlim = c(0,5),legend = rownames(t), space = 0.6,cex.names = 0.8,axis.lty = 1))

}

## [1] 0.22 0.54 0.86 1.18 1.50 1.82 2.14

## [1] 0.22 0.54 0.86 1.18 1.50 1.82 2.14

Contingency tables for barplots. Note that this makes more sense if done with only 2 samples, as only one statistic will be calculated per cluster.

clustercolumn = "seurat_clusters"
samplecolumn = c("sample_description")

all = table(sobject[[clustercolumn]][,1],sobject[[samplecolumn]][,1])
array = as.data.frame(cbind("Cluster" = 1, "Sample1" = 1, "Sample2" = 1, "p-value" = 1, "Cramer's V" =1))
for (i in 1:(ncol(all)-1)){
  for (k in (i+1):ncol(all)){
    partial = all[,c(i,k)]
    total = colSums(partial)
    clusters = rownames(partial)
    test = list()
    chi = chisq.test(partial)

    for (cluster in clusters){
      x = rbind(partial[cluster,],total-partial[cluster,])
      rownames(x) = c(cluster,"total")
      chi = chisq.test(x)
      f = sqrt(chi$statistic/sum(x))
      test[[cluster]] = chi
      test[[cluster]]$f = f
      
      pval = test[[cluster]]$p.value
      Cramer = sqrt(test[[cluster]]$statistic/sum(all[cluster,]))
      out = c(cluster,colnames(partial),pval,Cramer)
      #print out the p-value and cramer's v (closer to 1 = higher effect size) for each cluster
      array = rbind(array,out)
    }
  }
}
array = array[2:nrow(array),]
array

##    Cluster     Sample1     Sample2              p-value           Cramer's V
## 2        0  Control-P5  Control-P7    0.361085177576349   0.0109394808198687
## 3        1  Control-P5  Control-P7 7.64665654489035e-22    0.191938177732995
## 4        2  Control-P5  Control-P7 5.95614534513094e-29    0.315832817634304
## 5        3  Control-P5  Control-P7 0.000510828865315769    0.100148205635913
## 6        4  Control-P5  Control-P7    0.992014330250716  0.00059079605512425
## 7        5  Control-P5  Control-P7   0.0283761617600838    0.188661309213492
## 8        6  Control-P5  Control-P7    0.457057004447663   0.0968217050180083
## 9        0  Control-P5 Deprived-P5 1.34437556296484e-07   0.0631561954152618
## 10       1  Control-P5 Deprived-P5    0.240570556886315   0.0234524831429395
## 11       2  Control-P5 Deprived-P5   0.0061179336870449   0.0775382950903311
## 12       3  Control-P5 Deprived-P5    0.313212586343204    0.029064601013856
## 13       4  Control-P5 Deprived-P5 5.74325465511809e-61    0.972365568980841
## 14       5  Control-P5 Deprived-P5    0.539737740038465    0.052776681627256
## 15       6  Control-P5 Deprived-P5    0.384777394857507    0.113150763165275
## 16       0  Control-P5 Deprived-P7    0.909627625387819  0.00135959959488482
## 17       1  Control-P5 Deprived-P7 2.18972382631473e-26    0.212405620444991
## 18       2  Control-P5 Deprived-P7 5.61797419449398e-55     0.44170221339753
## 19       3  Control-P5 Deprived-P7    0.923147899499962  0.00278019739973707
## 20       4  Control-P5 Deprived-P7    0.966287651239758  0.00249480889152806
## 21       5  Control-P5 Deprived-P7 3.36736328288764e-06    0.399953288782341
## 22       6  Control-P5 Deprived-P7    0.254690956627095    0.148288976969146
## 23       0  Control-P7 Deprived-P5 5.46787971596413e-11   0.0785468668922356
## 24       1  Control-P7 Deprived-P5 2.58495939617642e-15    0.158062913908414
## 25       2  Control-P7 Deprived-P5 1.04460165868827e-14     0.21874223664928
## 26       3  Control-P7 Deprived-P5   0.0289975590096986   0.0629279919150176
## 27       4  Control-P7 Deprived-P5 7.72486774798787e-83     1.13813382192853
## 28       5  Control-P7 Deprived-P5    0.192595588304233    0.112139134298106
## 29       6  Control-P7 Deprived-P5   0.0769274145317663    0.230278022713877
## 30       0  Control-P7 Deprived-P7    0.297711436963763   0.0124733142024764
## 31       1  Control-P7 Deprived-P7   0.0605471021946722   0.0375057751916102
## 32       2  Control-P7 Deprived-P7 1.22108212629581e-11    0.191702166205112
## 33       3  Control-P7 Deprived-P7    0.001169424801631   0.0935549222150577
## 34       4  Control-P7 Deprived-P7    0.999999999999931 5.08903694414157e-15
## 35       5  Control-P7 Deprived-P7  0.00390611973253633    0.248356813635614
## 36       6  Control-P7 Deprived-P7     0.71897617142499   0.0468455443258497
## 37       0 Deprived-P5 Deprived-P7 4.77246258376725e-07   0.0603120282246738
## 38       1 Deprived-P5 Deprived-P7  1.9647586227062e-19    0.180160327635645
## 39       2 Deprived-P5 Deprived-P7 5.87760093287866e-38    0.364284993161216
## 40       3 Deprived-P5 Deprived-P7    0.393022579752167   0.0246161575174639
## 41       4 Deprived-P5 Deprived-P7 2.43986306811705e-57    0.942070014481652
## 42       5 Deprived-P5 Deprived-P7 0.000197431955903215    0.320362849330336
## 43       6 Deprived-P5 Deprived-P7   0.0368052047501156    0.271823810217491

write.csv(array,file = file.path(dir,"Spreadsheets",paste0(name,"_contingency.csv")),row.names = T,quote = F)